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Abstract 


In this paper, we consider the perturbation analysis for the periodic generalized coupled Sylvester 
(PGCS) equation. The normwise backward error for this equation is first obtained. Then, we present 
its normwise and componentwise perturbation bounds, from which the normwise and effective condition 
numbers are derived. Moreover, the mixed and componentwise condition numbers for the PGCS 
equation are also given. To estimate these condition numbers with high reliability, the probabilistic 
spectral norm estimator and the statistical condition estimation method are applied. The obtained 
results are illustrated by numerical examples. 
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1. Introduction 

In this paper, we consider the following matrix equation: 



( 1 . 1 ) 


where Ak, Ck £ R mxm , , D/~ £ R nxn , and Ek, Tfc £ R mxrl are ^he gj ven coefficient matrices, and 
Xk, Yfc £ R mx " are the unknown matrices satisfying X p+ \ = X\. Hereafter, R mxn denotes the set of 
m x n real matrices. 

The equation (1.1) is called the periodic generalized coupled Sylvester (PGCS) equation with pe¬ 
riod p (see e.g.,[3, 14]). It is easy to find that if p = 1, the PGCS equation reduces to the generalized 
coupled Sylvester (GCS) equation, which plays an important role in the linear control systems (see 
e.g., [5, 24]). One of the significant applications of this equation originates from computing the stable 
eigendecompositions of matrix pencils [6]. Some numerical methods were provided to compute the 
solution of the GCS equation (see e.g., [9, 18, 19]). Considering the specific structure of this equa¬ 
tion, Kagstrom [20] investigated its perturbation analysis, and derived the normwise backward error, 
normwise perturbation bounds, and normwise condition number. The derived results generalized the 
corresponding ones for the classic Sylvester equation given in [16]. Since the normwise condition num¬ 
ber cannot accurately reflect the influence of perturbations for some small entries in the data and 
ignores the structures of both input and output data with respect to scaling, Lin and Wei [26] pre¬ 
sented the mixed and componentwise condition numbers for the GCS equation. These two condition 
numbers were named by Gohberg and Koltracht [10]. The former measures the errors in output using 
norms and the input perturbations componentwise, and the latter measures both the errors in output 
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and the perturbations in input componentwise. To estimate the normwise, mixed and componentwise 
condition numbers for the GCS equation effectively, Diao et al. [7] applied the statistical condition 
estimation (SCE) method, which was proposed by Kenney and Laub in [21] and found applications in 
estimating the condition numbers of linear systems, least squares problem, eigenvalue problem, and 
matrix equations (see e.g., [7, 8, 13, 22, 23, 25]). Moreover, the authors also derived the effective 
condition numbers for the GCS equation and the classic Sylvester equation in [7], which can be much 
tighter than the normwise ones in [16, 20] in some cases. 

The PGCS equation also finds applications in many areas. For example, it can be used for structural 
analysis of periodic descriptor systems [4, 28]. Also, we will encounter this equation in computing 
periodic deflating subspaces associated with a specified set of eigenvalues [12]. So, some scholars 
considered the numerical methods for computing the solution of the PGCS equation, see e.g., [3, 14] 
and references therein. It was also shown in [12] that if 

A({(A fc ,G fe )}?)nA({( J B fe ,7? fe )}f) = 0, 

then the PGCS equation (1.1) has a unique solution. Here, A({(Gfc, Hk )}\) denotes the eigenvalue set 
of the periodic regular matrix pairs {{Gk, t?fc)}i . This condition is equivalent to the fact that the 
coefficient matrix of the matrix-vector form of (1.1) is nonsingular. The matrix-vector form is 

Wz = g, (1.2) 

where 
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and 

2 = vec ([Ad, Yi, • • • , X p , Y p \ ), g = vec ([Ed, Fi, ■ ■ ■ , E p , F p ]) . 

In the above expressions, X 0 Y denotes the Kronecker product [11], the operator ’vec’ stacks the 
columns of a matrix one underneath the other [11], / is the identity matrix of appropriate order, and 
K t stands for the transpose of the matrix K. 

For the similar motivations in [7, 8, 16, 20, 26], we investigate the perturbation analysis for the 
PGCS equation in this paper. After introducing the notation and preliminaries in Section 2, we 
present the normwise backward error for the PGCS equation in Section 3. In Section 4, the normwise 
and componentwise perturbation bounds for the PGCS equation are derived. A normwise condition 
number and the effective condition number are also given in this section. In Section 5, we provide 
the mixed and componentwise condition numbers for the PGCS equation. An algorithm based on 
the SCE method is proposed to estimate the mixed and componentwise condition numbers in Section 
6. To estimate the normwise and effective condition numbers, we consider an alternative method, 
that is, the probabilistic spectral norm estimator by Hochstenbacli [17], which provides a reliable 
estimation of the spectral norm. A corresponding algorithm is devised in Section 6. In addition, the 
numerical examples are also given in this section to illustrate the differences between the normwise, 
effective, mixed and componentwise condition numbers, and the efficiency of the statistical condition 
estimations, respectively. Finally, we present the conclusion of the whole paper. 

2. Notation and preliminaries 

For the matrix A = ( a,ij ) € W mxn , Al, ||A|| 2 , UAUoo, and ||A||^ stand for its Moore-Penrose inverse, 
spectral norm, max row norm, and Frobenius norm, respectively, |A| is the matrix with elements 
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\atj\, and ||^|| max is defined by ||A|| max = ||vec(A)|| 00 . For the vectors a = ,a p ] T £ R p and 

b = [bi, ■ ■ ■ , b p ] T £ R p , we define the entry-wise division between a and b by a/b = [ci, ■ • • , c p \ T with 

_ f 5 if 6* 7^ 0, 

1 \ a u if bi = 0. 


Following [29], the componentwise distance between a and b is defined by 


d(a , b ) = 


a — b 

— max < 

’ \a% ~ 6*| ' 

1 [ 

b 

*=!,••• 1 P 

OO 

l N J 



if K ^ o, 

if b io = 0. 


Note that when bi 0 0, d(a, b) gives the relative distance from a to b with respect to b , while the 
absolute distance for bi 0 =0. 

In order to define the mixed and componentwise condition numbers, we also need to define the 
set B°(a,e) = {x = [xi, ■ ■ ■ , x p ] T £ R p | |x* — a*| < e|a*|, * = 1, • • • ,p} with a = [ai, • • • , a p ] T £ R p 
and e > 0, and denote the domain of definition of a function F : R p —> R 9 by Dorn(F). Thus, the 
definitions of the mixed and componentwise condition numbers can be given as follows. 

Definition 2.1. [29] Let F : R p —> R 9 be a continuous map defined on an open set Dom(F) C R p 
such that 0 ^ Dom(F'). Let a £ Dom(F') ; such that F(a) 7^ 0. 

1. The mixed condition number of F at a is defined by 


m(F , a) = lim sup 

e_>0 x^a 

x6B°(a,e) 


||F(x)-F(a)|| 00 1 

IIWIL d{x,a)' 


2. The componentwise condition number of F at a is defined by 

d(F(x),F(a)) 


c{F , a) = lim sup 

x -£ a 

xeBO(a, e ) 


d(x, a) 


The Frechet derivative is essential in deriving the explicit expressions of condition numbers. Its 
definition is presented below. 

Definition 2.2. Let X and y be two Banach spaces, and a map F : U £ X —> y with U being an 
open set. Then F is said to be Frechet differentiable at a £ U, if there exists a bounded linear operator 
DF a : X —> y such that 

\\F(a + h)-F(a)^DF a (h)\\ 

h^o ||/i|| 

When the map F in Definition 2.1 is Frechet differentiable, the following lemma given in [29] reduces 
the computation burden of mixed and componentwise condition numbers. 


Lemma 2.1. Under the assumptions of Definition 2.1, when F is Frechet differentiable at a, we have 


m{F, a) 


c(F, a) 


I DF(a ) || a | 

IIWIloo 

I PF(a) || a I 
I F{a) | 


where DF(a ) is the Frechet derivative of F at a. 


( 2 . 1 ) 

( 2 . 2 ) 


To estimate the mixed and componentwise condition numbers, we need the SCE method which is 
ever mentioned in Section 1. In the following, we present a brief introduction on this method. 
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For a twice continuously differentiable function / : K. p —► R, by Taylor’s theorem, we get 


/(x + 5z) = f(x) + SXf(x) T z + 0{5 2 ), 


(2.3) 


where <5 is a small positive number, V/(x) = 


9f(x) df(x) 
dx± 5 8 x 2 ’ 


9f(x) 
5 dx v 


is the derivative of / at x, and 


2 £K p satisfies \\z \\2 = 1. From (2.3), the following inequality can be derived easily 

\f(x + 6z ) - f(x)\ « <5|V/(x) T z| < tf||V/(a:)|| 2 , 


which shows that the local sensitivity can be measured by a magnification factor S and the absolute 
condition number || V/(a?) || 2 - Based on the firm theoretical analysis given in [21], we have that if we 
choose a random vector z from U(S P - 1 ), the uniform distribution over unit sphere S p -1 in R p , then 
the following equality holds 

E(| X f{x) T z |) = Wp||V/(a;)|| 2 , (2.4) 

where E(-) is the expectation operator, and u> p is the Wallis factor with u)\ = 1, w 2 = 2/ 7 r, and 


UJ p — 


1- 3-5 —(p—2) 

2- 4-6—(p-1) ’ 

2 2-4-6---(p—2) 
7r 3-5-7—(p-1) ’ 


for p odd, 

when p > 2 . 

for p even, 


Owing to the equality (2.4) and the easy approximability of the Wallis factor (u p ~ ^(p -±) P reserves 

high accuracy), rj = \ V f(x) T z |/w p can be used as a condition estimator, and satisfies the following 
probability relationship 

p r ( l|V/(.r )|| 2 ^ < -y||V/(aj)|| 2 ) > 1 - — + 0{\), with 7 > 1 . 

7 7T7 7 Z 

According to [21], the accuracy of condition estimator can be enhanced by multiple samples. If we 
choose two samples z\, Z 2 G U(S P - 1 ), then the condition estimator given by 


? 7( 2 ) = — J\V f{x) T Zi\ 2 + |V/(x) t 2 2 | 2 

OJp v 

with Zi, Z 2 being obtained from Z\ and z 2 by orthonormalization meets the following probability 
relationship 


Pr ( l|V/(a ' )l12 < r?(2) < 7 ||V/(z)|| 2 ) 


1 - 


47 2 


In the similar manner, a general fc-sample SCE estimator can be defined [21]. 

In addition, in the following sections, we will apply the following equality frequently 


vec (AXC) = (C T <8> A)vec(X), 


(2.5) 


where A, X and C are matrices of appropriate orders such that the product AXC is well-defined. The 
equality (2.5) can be found in [11]. 


3. Normwise backward error 


Let Z = 


Xi ,Y 1} - 


,X P ,Y P 


denote an approximate solution to the PGCS equation (1.1). The 
normwise backward error of Z is defined by 

rj{Z) = min{e : (A k + A A k )X k - Y k {B k + A B k ) = E k + A E k , 

{■C k + A C k )X k+1 - Y k (D k + A D k ) = F k + A F k , k = !,-■■ ,p}, (3.1) 
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where AAk,ACk G 


ABfc, A Dk G 


l , and AE k , AF k G 


satisfy 


||Aj4fe||jr < eat, ||ABfc||i? < e/3fc, ||AU fe ||F < ey*,, ||ACfc||ir < eCfc, \\AD k \\F < et*,, ||AFfc||i? < edfc. 

(3.2) 

The tolerances a*,, p k , 7 *,, and d*, provide some freedom in how we measure the perturbations. 

Usually, 

CKfe = ||j4fc||F, Pk = ||-Bfc||^, 7fc = H-E^Hf, Cfc = IICfcllF) Tfc = ||Dfc||F, dk = 11-PfcllF- (3-3) 

In this case, the normwise backward error is called the relative normwise backward error with respect 
to Frobenius norm. 

The equation in (3.1) can be rewritten as 


AAkXk — Y k A B k — A E k — E k — {A k X k — YkB k ) — R k \, 
ACkXk+i — YkADk — A Fk = Fk — (CkXk +1 — Y k D k ) = Rk2, 


k = !,■■■ ,p, 


(3.4) 


where 1Z = [Rn,Ri 2 , , R P i, Rp'i] denotes the residual corresponding to the solution Z. Using the 
Kronecker product and (2.5), we can rewrite (3.4) as 


(A'J 0 J)vec(AA fc ) - (I 0 Y k )vec(AB k ) - vec(A E k ) = vec(R k i), 
{X k+1 0 I)vec(AC k ) - (10 Y k )vec(AD k ) - vec(AFfe) = vec(i? fc2 ), 


k = !,-■■ ,p- 


That is, 


Hu = r, 

where H = diag(7/i, • • • , H p ) with 

a k(X£ 01) ~p k (l0Y k ) - 7 k I 


(3.5) 


H k = 


0 0 0 
C k (X k+1 0l) -T k (l0Y k ) -5 k I 


and 


u = 


vec(AAi) T vec(AUi) T vec(AFd) T vec(AC'i) 7 ’ vec(ADi) T vec(AFi)’ ] 


OL l 


Pi 


7i 


Cl 


Tl 


<5i 


vec(AA p ) T vec(AB p ) T vec(AE p ) T vec(AC P ) T vec(AD p ) T vec(AU p ) 5 


T 


Op Pp 

r = vec ([i?n, i? 2 i, • • • ,i?pi,i? P 2 ]) ■ 


7p 


Cp 


It is easy to find that 17 is full row rank if 7 fc / 0 and 5 k 7 ^ 0 for k = 1, • • • ,p. In this case, (3.5) has a 
minimum Euclidean norm solution 

u = H'r. 

From the definition of normwise backward error, we have 


v{z) < 


Wi 


On the other hand, considering (3.2) 

IA B 


11 112 i|AAj|| F 

Nla = E-3- 

1=1 1 


+ 


• IlF , \\ AE i\\ F , IIACill 

j-r - 9 - ■ 


IA D 


i\\F 


IA Fi 


* 11 F 


Therefore, 






Ti 


< V(Z) < 


c 2 


iff 7 




< 6 pe. 


(3.6) 


Thus, we obtain both the upper and lower bounds of the normwise backward error r](Z) for the PGCS 
equation. 
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Remark 3.1. If the period p = 1, the bounds in (3.6) reduce to the corresponding ones for the GCS 

equation. The reduced lower bound is a little different from the one in [20] since the definitions of 

normwise backward error here and in [20] are a little different. Further, if C\ = 0, D\ = 0, and F\ = 0, 
we have the results for the classic Sylvester equation [16]. Note that \/6 should be replaced by a/ 3 in 
this case. 

4. Perturbation bounds 

Assume that the matrices A k , B k , E k , C k , D k ,F kl X k and Y k in (1.1) are perturbed as 

Ak —> Ak + A Ak, B k —> B k + AB k , E k — > Ek + A Ek, Xk Xk + XX k , 

Ck —t Ck + A Ck, Dk —> Dk + A Dk, Fk Fk + A Fk, Yk —¥ Yk + A Yk, 

where AA kl AC k £ M mxm , AB k ,AD k £ K nx ", AE k , AF k , AX k , AY k £ K mxn , and AX P+1 = AX ± . 
Then the perturbed PGCS equation (1.1) is 

f ( A k + A A k )(X k + AX k ) — ( Y k + AY k ){B k + AB k ) = E k + A E k , 

< k = l, ■■■ ,p. (4.1) 

{ (C k + AC k )(X k +i + AX fc+1 ) - (Y k + AY k )(Dk + AD k ) = F k + AF k , 

In the following, we regard AX k . AY k (k = 1, • ■ • ,p) as the unknown matrices of the matrix equation 
(4.1), and obtain the condition under which the equation (4.1) has the unique solution, and then the 
desired perturbation bounds. 

Considering (1.1), the equation (4.1) can be simplified as 

( A k AX k — AY k B k = AE k — (AA k X k — Y k AB k ) — (AA k AX k — AY k AB k ), 

\ CkAX k+1 - AYkDk = AFk - (A C k X k+ i - Y k AD k ) - (AC k AX k+ i - AY k AD k ), 

which, using the Kronecker product and (2.5), can be rewritten as 



vec(AXi) 


vec(A£i) 


vec(Xi) 


vec(AAi) 

w 

vec(AFi) 

_ 

vec(AFi) 

- AW 

vec(Fi) 

- AW 

vec(AFi) 


vec(AA p ) 


vec(A E p ) 


vec(Xp) 


vec(AX p ) 


vec(AYp) 


vec(A F p ) 


vec (Y p )_ 


vec(A Y p ) 


where AW is the same as W in (1.2) with A k ,B k , C kl and D k being replaced by AA kl AB k , AC kl and 
A D k , respectively. Let 

A^ = vec ([AXi, AYi, ■ • ■ , AX. p , AF p ]) , A g = vec ( [AEx , AFi, • • • , A E pi AF P ]) . 

Then we simplify (4.2) as 

WAz = Ag - AWz - AWAz. (4.3) 

Combining the first two terms in the right side of (4.3), we can rewrite (4.3) as 

WAz = -H lU - AWAz, (4.4) 

where H\ is the same as H in (3.5) except that X k and Y k in (3.5) are replaced by X k and Y k , 
respectively. Thus, 

A^ = - W-'AWAz. 

Define the operator equation of A^; as follows 

$(Az) = Ai = —W~ l H\u - W _1 AWA*. (4.5) 
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In the following, we use the Banach fixed point theorem (see, e.g., [24, Appendix D]) to derive the 
bound for A z. 

Let 

||W r_1 AW'|| 2 < 1, 

and denote the set f l as 


(4.6) 


12 = < s £ K. 2mrip : ||s|| 2 < 


Iw-'HM 


{ i-\\w-'aw\\ 2 ' 

which is closed and convex. Then, for any si,s 2 £ 12, we have 

||$(<»i)|| 2 < ||W r - 1 J ffi«|L H- ||W r - 1 AW r |L ||si|| 2 


W- l H lU \ 


„ 1 II II 1 „ \\ W ~ lH M\o 

< \\W~ 1 Htu\\ +\\w~ 1 aw\\ — -—— = — -—— 

- 11 1 1,2 + 11 1,2 1 - WW-iAWWz l-WW-'AWWz'- 


and 


||$(si) - 4>(s 2 )|| 2 < WW-'AWs! - W- 1 AWs 2 \L < WW-'AWL ||si - s 2 || 2 . 


Therefore, <&(•) maps the set 12 into itself and is contractive (see, e.g., [24, Appendix D[). According 
to the Banach fixed point theorem, we have that there is a unique solution A z to the equation (4.5) 
in the set 12 when (4.6) holds. As a result, 


||Az|| 2 = ||[AX,AF 1 ,... ,AX p ,AF p ]|| f < ^ 
\\W~ 1 Hi | 


W^Hml 


< 


i—l * 


HAS; 


\\w- i aw\\ 2 

IIAE.ld 


l|AC»r p 

C? 


IIADi 


l|AFi| 


1/2 


What’s more, if set 


1-||W-1AW|| 2 

lAAiHf IIABillf ||A£ 1 || f IIACiHf \\AD4f ||AFi||f 


(4.7) 


e = max 


cti /3i 71 Ci n 8 \ 

|AA p || f ||AB p ||f ||AL; p || f ||AC p || f IIAI? p ||f \\AF p \\ f 


Pp 7 p Cp 


then we have 


v / 6p||W / - 1 I7i|l e 

||[A.Y,Ar 1 ,--. ,a.x p ,ay p ]\\ f < i - » w _ lM If . 


(4.8) 


In summary, we have the following theorem. 


Theorem 4.1. Assume that the unperturbed and perturbed PGCS equations are given in (1.1) and 
(4.1), respectively. If the perturbations in (4.1) satisfy (4.6), then the perturbed PGCS equation (4.1) 
has a unique solution, and the normwise perturbation bounds (4.7) and (4.8) hold. 

Remark 4.1. From (4.4) or (4.8), by omitting the high-order terms, we can get the following first- 
order perturbation bound 


|| [AX, AFi, • • • , AX P , AY P ]|| F < ^ HIT" 1 #! \ 

The above bound is attainable to first-order in e. So, 

\\w- x h 1 \L 


k N l = 


||[x 1 ,y 1 ,--- ,x p ,Y r 


(4.9) 


(4.10) 


pi 4pJ||F 


can be regarded as the normwise condition number for the PGCS equation (1.1). It is a generalization 
of the ones for the GCS equation and the classic Sylvester equation given in [16, 20]. 


7 



























Remark 4.2. Using the equation (4.3), along the same line for deriving (4.7), we have the following 
bound under the condition (4.6), 


||[AX 1 ,AFi,--- ,AX p: AY plUF 


\\[Xi,Yi,--- ,X p ,Y f 


Pi II F 


< 


< 


\\ W 1 112 IA 3 II 2 I llTIA — 1 A TJA 

||[.Y 1 ,y 1 ,..,Y p ,y P H! F + IF AW 


1 - ||W-iAW|| 2 

||[AE 1 ,AF 1 ,-,AF p ,AF p ]|| f L , \\&W\\2 

111 E u F 1 .-,E p ,F p \\\ f l 'E+ lmr ^W) 


1 — ||W r - 1 AW r || s 


where k{W) = ||IU|| 2 ||IU 1 1 |2 and 


kE = 


||fU- 1 || 2 ||[£; 1 ,F 1 ,... ,E p ,F t 


P\\\F 


||[x 1 ,F 1 ,-.. ,x p ,v r 


PJIlF 


As done in [7], we can call kE the effective condition number for the PGCS equation (1.1). It can be 
much tighter than fc/vi if there are only perturbations on the right-hand side of the equation (1.1). 
The main reason is that kE only contains the information of [E\, T\, • • • ,E p , F p ], while fc/vi contains 
the information of all the coefficient matrices. 


Now we consider the componentwise perturbation bounds for the PGCS equation using the operator 
equation (4.5) and the generalized Banach fixed point theorem (see, e.g., [24, Appendix D]). 

Let 


radius (|W _1 AW|) < 1, 


(4.11) 


and define the set 3 as 


3= {ser np : |s| < (/- |IU- 1 AlU)|)- 1 |IU- 1 ffiu|} . 

It is easy to check that the set 5 is closed and convex, and for any si, s 2 £ 3, 

|$(si)| < |W - 1 J?iu| + |W _ 1 AW| |si| 

< \W~ l Hiu\ + |IU _ 1 AIU| (/- |IU - 1 AIU )|) -1 |W _ 1 ITiu| 

= (/-|IU" 1 AIU)|)" 1 |IU" 1 Riu| 

and 

|$(si) - $(s 2 )| < {W-'AWs! - W~ 1 AWs 2 \ < |W _ 1 AW| |si - s 2 |. 

Therefore, 4>(ti, •) maps the set 3 into itself and is generalized contractive (see, e.g., [24, Appendix D]). 
According to the generalized Banach fixed point theorem, we have that there is a unique solution Az 
to the equation (4.5) in the set 5 when (4.11) is satisfied. As a result, 

|As| = vec([|AA 1 ,AYi,--- ,AX P ,AY P \]) < (I — |IU - 1 AIU )|) -1 \W~ 1 Hiu\. (4.12) 

The above discussions imply the following theorem. 

Theorem 4.2. Assume that the unperturbed and perturbed PGCS equations are given in (1.1) and 
(4.1), respectively. If the perturbations in (4.1) fulfill (4.11), then the perturbed PGCS equation (4.1) 
has a unique solution, and the componentwise perturbation bound (4.12) holds. 

Remark 4.3. Form (4.12), we have the first-order componentwise perturbation bound 

vec ([|AA'i, AFi, • • • ,AA P ,AF P |]) < (4.13) 

Remark 4.4. When the period p = 1, the perturbation bounds obtained in this section reduce to the 
corresponding ones for the GCS equation, where the first-order normwise one is equivalent to the one 
in [ 20 ] in essence. 













5. Mixed and componentwise condition numbers 

In this section, using Lemma 2.1, we investigate the mixed and componentwise condition numbers 
for the PGCS equation 

We first rewrite (4.4), omitting the high-order terms, as follows, 

WAz k, -H 2 v, 

where H 2 and v are the same as H i and u in (4.4), respectively, except that all the tolerances 
otk)Pk,~tk,Ck,Tk and 5 k are replaced by 1. Thus, 

Az —W~ 1 H 2 v. (5.1) 


Define the map 4/ as 


SP :t-> z, 


where t = [vec(7i) T , vec(I?i) r , vec(I?i) T , vec(Ci) T , vec(Di) T , vec(Fi) T ,-• ■, vec(7 P ) T , vec (B P ) T , 
vec(E p ) T , vec (C p ) T , vec (D p ) T , vec {E p ) T ] T : and z is defined as in (1.2). Then from Definition 2.2 and 
(5.1), it follows that the Frechet derivative of at t is: 


M(f) = -W~ l H 2 . 


(5.2) 


Thus, combining Lemma 2.1 with (5.2), we have the following theorem which gives the expressions of 
the mixed and componentwise condition numbers of the PGCS equation (1.1). 

Theorem 5.1. With the above notation , the mixed and componentwise condition numbers of the PGCS 
equation ( 1 . 1 ) are given by 


m(4/,t) = 
cOM) = 


\W~'H 2 \ \t\\ 


\\[XuY lr .. ,A"p,r p ]|| max 


||[x 1 ,Yi,--- ,Ap,y p ]|| max ’ 


| w- 1 ^! |t| 


UJ 

vec([ X U Y U ---,X P ,Y P ]) 

OO 

vec([X 1 , Y), ■ ■ • ,X p ,Yp\) 


(5.3) 

(5.4) 


where 


uj = 


IF " 1 

Af ® I 

0 

|vec(Ai)| + 

w- 1 

I 0 Y\ 

0 

|vec(Ri)| + 

w- 1 

71 

oj 




0 




0 




O' 


+ 

w- 1 

0 

* 


|vec(Ci)| + 

IF ” 1 

J 0 Yi 


vec(Di)| + 

IF ” 1 

I 




0 




0 




0 



vec(Si)| 
vec(Fi)| 




0 




0 




'O' 


_l-1- 

IF” 1 

0 

k 1 a 


|vec(7p)| + 

IF” 1 

/0Fp 


|vec(Rp)| + 

IF” 1 

I 




0 




0 




0 



|vec(£’ p )| 


+ 


W 


-1 


X? 0 I 


|vec(Cp)| + 


W 


-1 


I < 8 > Y p 


|vec(£) p )| 


W 


-1 


|vec(£’ p )|. 


Proof. In view of Lemma 2.1, (5.2) and the definition of ||-|| , it is necessary only to show how to 

obtain the expression of ui. This can be done easily by using the expressions of Pl 2 and t. □ 


Remark 5.1. Note that 


= |i||l < ||IF” 1 |I |||R 2 | mil = IIif- 1 ! 


|A 1 ||X 1 | + |Yi||Bi| + |f? 1 |' 
IGiH^I + IFiHDil + IFrl 

|Ap||X p | + |l p ||Rp| + \E p \ 

|Cp||Xi| + IFpllDpI + |Fp|j 
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Here, the definition of ||-|| max 


is used. So, we have an upper bound for the mixed condition number 


m('F, t) < 


\W~ 


||[X 1 ,Fi,... ,X p ,F p ]|| max 


]Ai||Xi| + lYiH-Bil + |Ui| 

|Ci||X 2 | + |y 1 ||D 1 | + |F 1 | 

\A P \\X P \ + \Y P \\B P \ + \E P \ 
l\C P \\Xi\ + \Y p \\D p \ + \F p \j 


From the definition of the entry-wise division of vectors given in Section 2, we have 

(diagtvecaXi.yi,--- ,X p ,Y p ]))f \W~ 1 H 2 \ |f| 


3 

1 

to 

C-+- 


vec([ X^Yu---,X P ,Y P }) 

OO 


(5.5) 


Here, for a vector a = [cti,--- , a p ] T € R p , (diag(a))* denotes a diagonal matrix with the elements 
aj(* = 1, •• • ,p) of the following form 

t f ^ ^ ^ 0, 

’ ' 1, if a n = 0. 


Then 


cOM) < 


(diag (vec([A'i, hi, • • • , X p , Fp])))* W 1 


IIIWIIL- 


oo 

Thus, an upper bound for the componentwise condition number can be given by 


lyillBil + 1^1 

\C 1 \\X 2 \ + milDil + |Fi| 


c(T,t) < 


(diag (vec([A'i, Hi, • • • , X p , F p ])))* W 1 


OO 


(5.6) 


\ a p\\ x p\ + l^pll-®pl + \Ep\ 

LlCpIMM + \ y p\\ d p\ + I f pU max 

Remark 5.2. Using (5.2) and the definition of the normwise condition number given in [27], we can 
obtain an alternative normwise condition number for the PGCS equation (1.1): 


jy- 1 #. 


= 


J2 + IIAIIf + \\ e 4f + IICill_F + IIA|| f + II 

2 = 1 


1/2 


||[x 1 ,y 1 ,-- - ,x p ,y p ] 


which is a little larger than fcjvi in (4.10) if the tolerances a*,, Pk, 7fc, Cfci T fc an( l in (4-10) are chosen 
as in (3.3). In addition, if the period p = 1, the above condition number reduces to the corresponding 
one for the GCS equation [26]. 


6. Numerical experiments 

In this part, our attention mainly focuses on the comparison and estimation of the condition 
numbers derived in the above sections. 

We first provide an example to compare the normwise, effective, mixed and componentwise condi¬ 
tion numbers. The example is taken from [3] with some modifications. 
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Example 6.1. For the PGCS equation (1.1), let the period p = 3, and the coefficient matrices be 



"1 

0 

01 


A 

0.3 

8' 


'0.1 

0.03 

9 ' 


'1 

0 

12' 

9 

Ai = 

0 

1 

10 

5 ^2 — 

0 

1 

10 

5 ^3 — 

0 

0.1 

0.9 

= 


0 

0 

1 


0 

0 

1 


0 

0 

0.1 


Z 


1 

'1 21 ' 


'1 1' 


'0 1 ' 


'2 O' 

,b 3 = 

0 10“* 

,£l = 

0 1 

0 10 

,E 2 = 

2 1 

5 8 

,e 3 = 

3 1 
0 2 



'0.1 

10 1 . 5 ' 


'l.l 3 

8 


'l 0.5 0.9' 


1 n 

c\ = 

1 

10 0.1 

,c 2 = 

0.2 5 

0.1 

,c 3 = 

1 0.1 0.9 

■Di 

= 

1 u 

1 9 


2 

0.3 0.1 


1 0.01 0.01 


1 

2 0.15 


1 z 


9 O 


’1 1 


' 1 O' 


'0 l' 


'2 O' 


d 2 = 

z y 

2 1 

,d 3 = 

1 1 

3 10~ T 

,ir = 

0.1 1 

,f 2 = 

2 1 

,f 3 = 

3 1 

5 






2 0 


5 8 


2 5 



with t, t £ {1, 3,5}. 


Upon some computations, the numerical results are exhibited in Table I. 


Table I: Comparison of condition numbers 



r = 1 

t = 1 

r = 1 

t = 3 

r = 1 

t = 5 

r = 3 
t = 3 

r = 3 
t = 5 

r = 5 
t = 5 

kN 1 

564.1934 

1.4085e+003 

1.3455e+003 

1.4065e+003 

1.3438e+003 

1.3438e+003 

kN2 

2.3429e+004 

5.8489e+004 

5.5874e+004 

5.8407e+004 

5.5803e+004 

5.5803e+004 

kE 

263.9046 

182.1415 

181.5541 

182.1423 

181.5566 

181.5567 

m{^,t) 

52.9059 

18.1312 

16.1057 

18.1240 

16.1058 

16.1058 

cOM) 

1.3318e+003 

260.1651 

269.9788 

120.0864 

119.9581 

119.9582 


From Table I, one can easily find that the effective, mixed and componentwise condition numbers 
behave well in most cases, while the normwise condition numbers k^ri and kpj 2 may highly overestimate 
the condition of the PGCS equation. Here, it should be pointed out that c(*P,t) may be very large if 
there are very small elements in the solution. This may be the reason why c(\I/, t) is so large for r = 1 
and t = 1. In this case, some distinction should be made to cope with this extremal case. We suggest 
the projection method proposed by Arioli et al. [1], and Cao and Petzold [2], but we will not go that 
far in this paper. 

In the following, we will devise two algorithms based on the probabilistic spectral norm estimator 
and the SCE method to estimate the normwise, effective, mixed and componentwise condition numbers. 
The former will be called the PCE method for short. 


Algorithm 1 PCE for the normwise and effective condition numbers 

1. Generate a starting vector vo from U{S q - 1 ) with q = 2 p(m 2 + n 2 + mn). 

2. Compute the guaranteed lower bound a and the probabilistic upper bound /3 of ||1F _1 Hi|| 2 (||VU —1 1| 2 ) by 
probabilistic spectral norm estimator (||W~ 1 Hi|| 2 ^ (3 (11VU^ 1 11 2 ^ /3) will hold with a given probability 
1 — e, where e is a user-chosen parameter). 

3. Compute the normwise and effective condition number by 


kpceNl — 


a + 0 


2||[A'i,Yi,-- - ,X P ,Y P ] 



® + A || {Ei, F \, ■ • ■ , E p , F p ) Hj, \ 

2 IKAhAV-- ,x p ,y p )\\ f J ■ 


The main part of Algorithm 1 is to estimate \\W 1 Hi|| 9 (||VF 1 1| 2 ) by probabilistic spectral norm 
estimator. A detailed analysis of the estimator was given in [17] by Hochstenbach. The author showed 
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Algorithm 2 SCE for the mixed and componentwise condition numbers 


1. Generate the random matrices (Jin, Ln, Mu, Su, Nu, Qn, ■ • ■ , Rpi, L v i, M p i, S p i, N p i, Q p i), •••, 
(Ri s , Lis, Mis, Sis, Ni s ,Qis, ■ ■ ■ , R ps , L ps , M ps , S P s, Np S ,Qp S ), where Rkj,Skj G R mxm , Lkj,Nkj G 
R" x ” and Mkj,Qkj G R mxn with k = 1, • • • ,p, j = 1, • • • , s, and all entries being in the standard 
normal distribution Xf(0, 1). Orthonormalize the matrix 


’vec(J?n) 

vec(Ln) 


vec(-Ri s )" 

vec(Lis) 


Lvec(Qpi) ••• vec(Q ps )J 

to get an orthonormal matrix [pi, • • • ,p s ], Then, convert pj into the matrix form 
(Rij , Lij , Mij , Sij , Nij. Qij, ■ - • , R p j , L p j , M p j , Spj i Npj , Q P j). 

2. Set q = 2 p(m 2 + n 2 + mn), get the approximates of u q and co s , and let 

(Rij , Lij , Mij , Sij , Nij , Qij , • • • , -Rpj , Lpj, Mpj , S'pj, TVpj 

i Qpi) 

= ( Rij, Lij, Mij, Sij, Nij, Qij, • ■ • , Rpj, Lpj, M P j, S P j, N P j ,Q P j ) 
o (Ai, Bi, Ei, Ci, Di, Fi ■ ■ ■ , A p , B p , E p , C p , D p , F p ). 

Here, the symbol o denotes the Hadamard product. 

3. For j = 1, • • • , s, solve the following PGCS equation 

f AkXkj F kjBk = M k j (RkjXk F kLkj ), 


CkN(kj-ljj F kjDk — Qkj (Skj Xk j -1 \kNkj), 


k =!,-■■ ,p, 


and compute the absolute condition vector 


LUs \ - 

KabS= ^\h 


where Uj = vec ([Xij, Yij, ■ ■ ■ ,Xpj,Y p j\). Here, the operations of taking square root and power are 
componentwise. 

4. Compute the estimations of the mixed and componentwise condition numbers by 


nw('M) = 


\\[Xi,Yi,--- ,X P ,Y P ] 


Csce(®,t) = 


vec([Xi, Yi, • • • ,X p ,Yp})\ 


Note: For the sake of convenience, we write {Ai,Bi,Ei,Ci,Di,Fi - ■ ■ , A p , B p , E P ,C P , D p , F p ) as a matrix 
though the matrices in the parenthesis do not have same orders. 





that ||^ _1 -ffi|| 2 (ll^ -1 ^) can be contained in a small interval [a,p\ with high probability. Here 
P/a ^ 1 + 5, where <5 is another user-chosen parameter. In our computation, we take e = 0.001 
and S = 0.01 . Thus, ||W _1 Iii|| 9 < P (||1T _1 ||2 ^ P) holds with a probability at least 99.9% and 
P/a < 1.01. Hence, we take (a + P )/2 as the estimation of ||lT _1 .ffi|| 2 (|| 1 1| 2 )■ 

For Algorithm 2, we would like to choose s = 3 in numerical experiments. This means that 
m sce (T, t ) and c sce ('k, t) fall into the intervals [0.2 x m(\ &, t),5x m(\H, f)] and [0.2 x c('k, t), 5 x c( v k, f)] 
with the probability 1 — 3/^3 ~ 0.9913, respectively, if 7 = 5. 

Now we present a specific example to investigate the efficiency of these two algorithms in estimating 
the condition numbers. 

Example 6.2. For the PGCS equation (1.1), let p = 3, m = 5, and n = 4, and generate the coefficient 
matrices as follows: Ak,Ck G randn(m), Bk,Dk G randn(n), and Ej~,Fk G randn(m,n). Here, 
the Matlab functions are used. Since the orders of the coefficient matrices are not so large, we get 
the solution by solving the linear equation (1.2). The computed solution z satisfies the inequality 
111 W 1 11 111 00 /11 ^ 11 00 < 10 -8 [15, p.131] and is treated as the exact solution. We test 1000 PGCS 
equations, and define the ratios of the estimated condition numbers and the exact ones as follows 

kpceNl kpceE ^sce(^j^) C sce (fIt, t) 

kjvi ke m(w,t) c(w,t) 

Upon computation, we have the numerical results of these ratios and their means and variances: 
E(rjvi) = 1.0003, V(r-jvi) = 5.7960e- 007, E(r s ) = 1.0004, V(r E ) = 8.0694e- 007, E(r m ) = 1.8313, 
V(r m ) = 2.4788, E(r c ) = 2.4269, V(r c ) = 7.1857. The numerical results are plotted in Figure 1. 


The efficiency of estimating k N1 The efficiency of estimating k £ 



The efficiency of estimating m('F,t) The efficiency of estimating cC'P.t) 




Figure 1: Efficiency of condition estimators 


From Figure 1 and the results on means and variances, we can find that both the PCE method and 
the SCE method can give reliable estimations of the normwise, effective, mixed and componentwise 
condition numbers, respectively. 

Remark 6.1. In Example 6.2, we get the solution to the PGCS equation (1.1) by solving the linear 
system (1.2). The cost will be very expensive when the orders of the coefficient matrices in the PGCS 


13 

















































































































































equation (1.1) are large. In this case, other iterative methods need to be consulted; see [3, 14] and 
references therein. 

7. Conclusion 

In this paper, we investigated the perturbation analysis of the PGCS equation. The normwise back¬ 
ward error for this equation is first given. Then, by Banach fixed point theorem, we derive its rigorous 
normwise and componentwise perturbation bounds, from which the first-order perturbation bounds, 
and the normwise and effective condition numbers are obtained. Moreover, the explicit expressions of 
the mixed and componentwise condition numbers and their upper bounds for the PGCS equation are 
also given. A simple example is provided to illustrate the differences among these condition numbers. 
To estimate these condition numbers, the probabilistic spectral norm estimator and the SCE method 
are introduced and two algorithms are devised. From the numerical experiments, we find that both 
the PCE method and the SCE method perform efficiently in estimating the normwise, effective, mixed 
and componentwise condition numbers, respectively. 
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